rm(list=ls(all=TRUE)) 


# Install and load packages
packages <- c("foreign", "Hmisc", "dplyr", "knitr", "kableExtra", "stargazer", "MASS", "lme4", "haven", "radiant", "psych",
              "magrittr", "lmtest", "multiwayvcov", "graphics", "xtable", "clubSandwich", "sjPlot", "sjmisc", "ordinal", "effects", 
              "lattice", "margins", "texreg", "foreach", "ggpubr", "sandwich", "pglm", "plm", "magick",  "list", "ggeffects",
              "dotwhisker", "broom", "oglmx", "optimx", "magick", "mlogit")

package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)


##########################################
##########################################
######## Mexico 2012 panel study  ######
##########################################
##########################################
mexico.data <- read.csv("Mexico 2012 panel study.csv", header=TRUE, stringsAsFactors=FALSE, row.names=NULL)

mexico.data$female_w1 <- mexico.data$p1
mexico.data$female_w1[mexico.data$female_w1==-1]<- NA
mexico.data$female_w1 <- mexico.data$female_w1 - 1

mexico.data$female_w2 <- mexico.data$w2_p1
mexico.data$female_w2[mexico.data$female_w2==-1]<- NA
mexico.data$female_w2 <- mexico.data$female_w2 - 1

mexico.data$panel <-0
mexico.data$panel[mexico.data$p1>0 & mexico.data$w2_p1>0 ] <- 1

mexico.data$refresh <- 0
mexico.data$refresh[mexico.data$p1==-1 & mexico.data$w2_p1>0] <- 1

mexico.data$family_benefit_a_w1 <-mexico.data$p38a
mexico.data$family_benefit_a_w1[mexico.data$family_benefit_a_w1== -1]<-NA
mexico.data$family_benefit_a_w1[mexico.data$family_benefit_a_w1== 3]<-NA
mexico.data$family_benefit_a_w1[mexico.data$family_benefit_a_w1== 9]<-NA
mexico.data$family_benefit_a_w1 <-2-mexico.data$family_benefit_a_w1

mexico.data$family_benefit_b_w1 <-mexico.data$p38b
mexico.data$family_benefit_b_w1[mexico.data$family_benefit_b_w1== -1]<-NA
mexico.data$family_benefit_b_w1[mexico.data$family_benefit_b_w1== 3]<-NA
mexico.data$family_benefit_b_w1[mexico.data$family_benefit_b_w1== 9]<-NA
mexico.data$family_benefit_b_w1 <-2-mexico.data$family_benefit_b_w1

mexico.data$family_benefit_c_w1 <-mexico.data$p38c
mexico.data$family_benefit_c_w1[mexico.data$family_benefit_c_w1== -1]<-NA
mexico.data$family_benefit_c_w1[mexico.data$family_benefit_c_w1== 3]<-NA
mexico.data$family_benefit_c_w1[mexico.data$family_benefit_c_w1== 9]<-NA
mexico.data$family_benefit_c_w1 <-2-mexico.data$family_benefit_c_w1

mexico.data$family_benefit_aggregate_w1<-rowSums(mexico.data[, c(which(colnames(mexico.data)=="family_benefit_a_w1" ), 
                                                                 which(colnames(mexico.data)=="family_benefit_b_w1" ), 
                                                                 which(colnames(mexico.data)=="family_benefit_c_w1" ))], na.rm = TRUE)
mexico.data$family_benefit_aggregate_w1 [ mexico.data$refresh ==1 ] <- NA

mexico.data$family_benefit_a_w2 <-mexico.data$w2_P39A
mexico.data$family_benefit_a_w2[mexico.data$family_benefit_a_w2== -1]<-NA
mexico.data$family_benefit_a_w2[mexico.data$family_benefit_a_w2== 3]<-NA
mexico.data$family_benefit_a_w2[mexico.data$family_benefit_a_w2== 9]<-NA
mexico.data$family_benefit_a_w2 <-2-mexico.data$family_benefit_a_w2

mexico.data$family_benefit_b_w2 <-mexico.data$w2_P39B
mexico.data$family_benefit_b_w2[mexico.data$family_benefit_b_w2== -1]<-NA
mexico.data$family_benefit_b_w2[mexico.data$family_benefit_b_w2== 3]<-NA
mexico.data$family_benefit_b_w2[mexico.data$family_benefit_b_w2== 9]<-NA
mexico.data$family_benefit_b_w2 <-2-mexico.data$family_benefit_b_w2

mexico.data$family_benefit_c_w2 <-mexico.data$w2_P39C
mexico.data$family_benefit_c_w2[mexico.data$family_benefit_c_w2== -1]<-NA
mexico.data$family_benefit_c_w2[mexico.data$family_benefit_c_w2== 3]<-NA
mexico.data$family_benefit_c_w2[mexico.data$family_benefit_c_w2== 9]<-NA
mexico.data$family_benefit_c_w2 <-2-mexico.data$family_benefit_c_w2

mexico.data$family_benefit_aggregate_w2<-rowSums(mexico.data[, c(which(colnames(mexico.data)=="family_benefit_a_w2" ), 
                                                                 which(colnames(mexico.data)=="family_benefit_b_w2" ), 
                                                                 which(colnames(mexico.data)=="family_benefit_c_w2" ))], na.rm = TRUE)
mexico.data$family_benefit_aggregate_w2 [ mexico.data$panel == 0 & mexico.data$refresh==0 ] <- NA

mexico.data$program_withdrawal_threat_w1 <- NA
mexico.data$program_withdrawal_threat_w1[mexico.data$family_benefit_aggregate_w1 >=0] <- 0
mexico.data$program_withdrawal_threat_w1[mexico.data$p39 ==1] <- 1
mexico.data$program_withdrawal_threat_w1[mexico.data$p39 == 3] <- NA
mexico.data$program_withdrawal_threat_w1[mexico.data$p39 == 9] <- NA

mexico.data$program_withdrawal_threat_w2 <- NA
mexico.data$program_withdrawal_threat_w2[mexico.data$family_benefit_aggregate_w2 >=0] <- 0
mexico.data$program_withdrawal_threat_w2[mexico.data$w2_P40 ==1] <- 1
mexico.data$program_withdrawal_threat_w2[mexico.data$w2_P40 == 3] <- NA
mexico.data$program_withdrawal_threat_w2[mexico.data$w2_P40 == 9] <- NA

mexico.data$personally_offered_w1 <- mexico.data$p40
mexico.data$personally_offered_w1[mexico.data$personally_offered_w1==-1] <-NA
mexico.data$personally_offered_w1[mexico.data$personally_offered_w1==3] <-NA
mexico.data$personally_offered_w1[mexico.data$personally_offered_w1==9] <-NA

mexico.data$personally_offered_w2 <- mexico.data$w2_P41
mexico.data$personally_offered_w2[mexico.data$personally_offered_w2==-1] <-NA
mexico.data$personally_offered_w2[mexico.data$personally_offered_w2==3] <-NA
mexico.data$personally_offered_w2[mexico.data$personally_offered_w2==9] <-NA

mexico.data$doubt_ballot_secrecy_w1 <- mexico.data$p36c
mexico.data$doubt_ballot_secrecy_w1[mexico.data$doubt_ballot_secrecy_w1 == -1] <- NA
mexico.data$doubt_ballot_secrecy_w1[mexico.data$doubt_ballot_secrecy_w1 == 5] <- NA
mexico.data$doubt_ballot_secrecy_w1[mexico.data$doubt_ballot_secrecy_w1 == 9] <- NA
mexico.data$doubt_ballot_secrecy_w1 <- mexico.data$doubt_ballot_secrecy_w1 -1

mexico.data$doubt_ballot_secrecy_w2 <- mexico.data$w2_P36C
mexico.data$doubt_ballot_secrecy_w2[mexico.data$doubt_ballot_secrecy_w2 == -1] <- NA
mexico.data$doubt_ballot_secrecy_w2[mexico.data$doubt_ballot_secrecy_w2 == 5] <- NA
mexico.data$doubt_ballot_secrecy_w2[mexico.data$doubt_ballot_secrecy_w2 == 9] <- NA
mexico.data$doubt_ballot_secrecy_w2 <- mexico.data$doubt_ballot_secrecy_w2 -1

mexico.data$age_w1 <- mexico.data$p2
mexico.data$age_w1[mexico.data$age_w1 ==-1] <-NA

mexico.data$age_w2 <- mexico.data$w2_p2
mexico.data$age_w2[mexico.data$age_w2 ==-1] <-NA

mexico.data$political_interest_w1 <- mexico.data$p5
mexico.data$political_interest_w1[mexico.data$political_interest_w1== -1] <- NA
mexico.data$political_interest_w1[mexico.data$political_interest_w1== 9] <- NA
mexico.data$political_interest_w1 <- 4 - mexico.data$political_interest_w1

mexico.data$political_interest_w2 <- mexico.data$w2_p5
mexico.data$political_interest_w2[mexico.data$political_interest_w2== -1] <- NA
mexico.data$political_interest_w2[mexico.data$political_interest_w2== 9] <- NA
mexico.data$political_interest_w2 <- 4 - mexico.data$political_interest_w2

mexico.data$state_w1 <- mexico.data$px6
mexico.data$state_w1[mexico.data$px6 == -1 ] <- NA

mexico.data$state_w2 <- mexico.data$w2_px6
mexico.data$state_w2[mexico.data$w2_px6 == -1 ] <- NA

mexico.data$urban_w1 <- mexico.data$px8
mexico.data$urban_w1[mexico.data$urban_w1== -1] <- NA
mexico.data$urban_w1 <- 2- mexico.data$urban_w1

mexico.data$urban_w2 <- mexico.data$w2_px8
mexico.data$urban_w2[mexico.data$urban_w2== -1] <- NA
mexico.data$urban_w2[mexico.data$urban_w2== 3] <- 1  
mexico.data$urban_w2 <- 2- mexico.data$urban_w2

mexico.data$religion_w1 <- mexico.data$p41
mexico.data$religion_w1[mexico.data$religion_w1==-1] <- NA
mexico.data$religion_w1[mexico.data$religion_w1== 9] <- NA

mexico.data$catholic_w1 <- NA
mexico.data$catholic_w1[mexico.data$religion_w1==1] <- 1
mexico.data$catholic_w1[mexico.data$religion_w1>1 & mexico.data$religion_w1 <=5] <- 0

mexico.data$other_religion_w1 <- NA
mexico.data$other_religion_w1[mexico.data$religion_w1>=2 & mexico.data$religion_w1<=4] <- 1
mexico.data$other_religion_w1[mexico.data$religion_w1==1 ] <- 0
mexico.data$other_religion_w1[mexico.data$religion_w1==5 ] <- 0

mexico.data$religion_w2 <- mexico.data$w2_P43
mexico.data$religion_w2[mexico.data$religion_w2==-1] <- NA
mexico.data$religion_w2[mexico.data$religion_w2== 9] <- NA

mexico.data$catholic_w2 <- NA
mexico.data$catholic_w2[mexico.data$religion_w2==1] <- 1
mexico.data$catholic_w2[mexico.data$religion_w2>1 & mexico.data$religion_w2 <=5] <- 0

mexico.data$other_religion_w2 <- NA
mexico.data$other_religion_w2[mexico.data$religion_w2>=2 & mexico.data$religion_w2<=4] <- 1
mexico.data$other_religion_w2[mexico.data$religion_w2==1 ] <- 0
mexico.data$other_religion_w2[mexico.data$religion_w2==5 ] <- 0

mexico.data$education_w1 <- mexico.data$p44
mexico.data$education_w1[ mexico.data$p44 == -1] <- NA
mexico.data$education_w1[ mexico.data$p44 == 10] <- NA

mexico.data$education2_w1 <- NA
mexico.data$education2_w1[mexico.data$education_w1==1 |mexico.data$education_w1==2 ] <-0
mexico.data$education2_w1[mexico.data$education_w1==3 |mexico.data$education_w1==4 ] <-1
mexico.data$education2_w1[mexico.data$education_w1>=5 & mexico.data$education_w1<=8 ] <-2
mexico.data$education2_w1[mexico.data$education_w1==9 ] <-3

mexico.data$education_w2 <- mexico.data$w2_P46
mexico.data$education_w2[ mexico.data$w2_P46 == -1] <- NA
mexico.data$education_w2[ mexico.data$w2_P46 == 10] <- NA

mexico.data$education2_w2 <- NA
mexico.data$education2_w2[mexico.data$education_w2==1 |mexico.data$education_w2==2 ] <-0
mexico.data$education2_w2[mexico.data$education_w2==3 |mexico.data$education_w2==4 ] <-1
mexico.data$education2_w2[mexico.data$education_w2>=5 & mexico.data$education_w2<=8 ] <-2
mexico.data$education2_w2[mexico.data$education_w2==9 ] <-3

mexico.data$economic_hardship_w1 <- mexico.data$p49
mexico.data$economic_hardship_w1[mexico.data$p49== -1] <- NA
mexico.data$economic_hardship_w1[mexico.data$p49== 9] <- NA
mexico.data$economic_hardship_w1 <- mexico.data$economic_hardship_w1 - 1

mexico.data$economic_hardship_w2 <- mexico.data$w2_P51
mexico.data$economic_hardship_w2[mexico.data$w2_P51== -1] <- NA
mexico.data$economic_hardship_w2[mexico.data$w2_P51== 9] <- NA
mexico.data$economic_hardship_w2 <- mexico.data$economic_hardship_w2 -1

mexico.data$partisanship_w1 <- NA
mexico.data$partisanship_w1[mexico.data$p22>=1 & mexico.data$p22 <=7] <-1
mexico.data$partisanship_w1[mexico.data$p23>=1 & mexico.data$p23 <=4 ] <-1
mexico.data$partisanship_w1[mexico.data$p23==0 ] <-0

mexico.data$partisanship_w2 <- NA
mexico.data$partisanship_w2[mexico.data$w2_P21>=1 & mexico.data$w2_P21 <=7] <-1
mexico.data$partisanship_w2[mexico.data$w2_P22>=1 & mexico.data$w2_P22 <=4 ]<-1
mexico.data$partisanship_w2[mexico.data$w2_P22==0 ] <-0

mexico.data$duration_community_residence_w1 <- mexico.data$p50
mexico.data$duration_community_residence_w1[mexico.data$p50 == -1] <- NA
mexico.data$duration_community_residence_w1[mexico.data$p50 == 99] <- NA 

mexico.data$duration_community_residence_w2 <- mexico.data$w2_P52
mexico.data$duration_community_residence_w2[mexico.data$w2_P52 == -1] <- NA
mexico.data$duration_community_residence_w2[mexico.data$w2_P52 == 99] <- NA 

# to combine
mexico.data.w1.sub<-mexico.data[c("personally_offered_w1", 
                                  "doubt_ballot_secrecy_w1",
                                  "family_benefit_aggregate_w1",
                                  "female_w1", 
                                  "age_w1", 
                                  "economic_hardship_w1",
                                  "political_interest_w1",
                                  "urban_w1",
                                  "panel", 
                                  "refresh",
                                  "religion_w1", 
                                  "catholic_w1",
                                  "other_religion_w1", 
                                  "education2_w1", 
                                  "state_w1", 
                                  "partisanship_w1", 
                                  "duration_community_residence_w1", 
                                  "program_withdrawal_threat_w1"
)]

for ( col in 1:ncol(mexico.data.w1.sub)){
  colnames(mexico.data.w1.sub)[col] <-  sub("_w1", "", colnames(mexico.data.w1.sub)[col])
}
mexico.data.w1.sub$wave <- 1

mexico.data.w2.sub<-mexico.data[c("personally_offered_w2", 
                                  "doubt_ballot_secrecy_w2",
                                  "family_benefit_aggregate_w2",
                                  "female_w2", 
                                  "age_w2", 
                                  "economic_hardship_w2",
                                  "political_interest_w2",
                                  "urban_w2",
                                  "panel", 
                                  "refresh",
                                  "religion_w2", 
                                  "catholic_w2",
                                  "other_religion_w2", 
                                  "education2_w2", 
                                  "state_w2", 
                                  "partisanship_w2", 
                                  "duration_community_residence_w2", 
                                  "program_withdrawal_threat_w2"
)]

for ( col in 1:ncol(mexico.data.w2.sub)){
  colnames(mexico.data.w2.sub)[col] <-  sub("_w2", "", colnames(mexico.data.w2.sub)[col])
}
mexico.data.w2.sub$wave <- 2

mexico.data.combined <- rbind(mexico.data.w1.sub, mexico.data.w2.sub)

###################
## Appendix A.2 ##
###################
temp.mexico.panel <- mexico.data[c("doubt_ballot_secrecy_w1",
                                   "personally_offered_w1", 
                                   "age_w1", 
                                   "female_w1", 
                                   "education2_w1", 
                                   "urban_w1",
                                   "economic_hardship_w1",
                                   "family_benefit_aggregate_w1",
                                   "program_withdrawal_threat_w1",
                                   "duration_community_residence_w1", 
                                   "catholic_w1",
                                   "other_religion_w1", 
                                   "partisanship_w1",
                                   "political_interest_w1",
                                   "doubt_ballot_secrecy_w2",
                                   "personally_offered_w2", 
                                   "age_w2", 
                                   "female_w2", 
                                   "education2_w2", 
                                   "urban_w2",
                                   "economic_hardship_w2",
                                   "family_benefit_aggregate_w2",
                                   "program_withdrawal_threat_w2",
                                   "duration_community_residence_w2", 
                                   "catholic_w2",
                                   "other_religion_w2", 
                                   "partisanship_w2",
                                   "political_interest_w2")]
sink("table_Appendix_A_2.tex")
stargazer(temp.mexico.panel, 
          type = "latex",
          summary.stat = c("n", "mean", "sd", "min", "max"),
          covariate.labels = c("Doubt about ballot secrecy", "Private transfer offered", 
                               "Age", "Female", "Education","Urban", "Economic hardship", "Num. of family benefit sources", 
                               "Threat of family benefit withdrawal", "Duration of community residence",
                               "Catholic","Other religion", "Partisanship", "Political interest",
                               "Doubt about ballot secrecy", "Private transfer offered", 
                               "Age", "Female", "Education","Urban", "Economic hardship", "Num. of family benefit sources", 
                               "Threat of family benefit withdrawal", "Duration of community residence",
                               "Catholic","Other religion", "Partisanship", "Political interest"
          ),                     
          title="Descriptive statistics of Mexico 2012 panel data", 
          digits=3 
)
sink()


###################
## Appendix A.3 ##
###################
mexico.panel.1<- polr(as.factor(doubt_ballot_secrecy) ~ personally_offered + 
                        wave +
                        panel + 
                        refresh +
                      as.factor(state)
                      , data = mexico.data.combined, Hess = "TRUE")

mexico.panel.2<- polr(as.factor(doubt_ballot_secrecy) ~ personally_offered + 
                        wave +
                        panel + 
                        refresh +
                        age +
                        female +
                        education2+
                        urban +
                        economic_hardship +
                        family_benefit_aggregate +
                        program_withdrawal_threat +
                        duration_community_residence +
                        catholic +
                        other_religion +
                        as.factor(state)
                      , data = mexico.data.combined, Hess = "TRUE")

mexico.panel.3<- polr(as.factor(doubt_ballot_secrecy) ~ personally_offered + 
                        wave +
                        panel + 
                        refresh +
                        age +
                        female +
                        education2+
                        urban +
                        economic_hardship +
                        family_benefit_aggregate +
                        program_withdrawal_threat +
                        duration_community_residence +
                        catholic +
                        other_religion +
                        partisanship +
                        political_interest +
                        as.factor(state)
                      , data = mexico.data.combined, Hess = "TRUE")

sink("table_Appendix_A_3.tex")
stargazer(mexico.panel.1, mexico.panel.2, mexico.panel.3,
          type = "latex",
          title="Regression results using Mexico 2012 Panel data", 
          omit = "state",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("State dummies", "Included", "Included", "Included")),
          covariate.labels = c("Private transfer offered", "Wave", "Panel respondent",  "New respondent", 
                               "Age", "Female", "Education","Urban", "Economic hardship", "Num. of family benefit programs", 
                               "Threat of family benefit withdrawal", "Duration of community residence",
                               "Catholic","Other religion", "Partisanship", "Political interest"), 
          notes.append = FALSE, 
          notes.align = "l")
sink()

##############################
###  Figure 1. Upper panel ###
##############################
predict.mexico.1 <- ggpredict(mexico.panel.3, terms=c("personally_offered"))
predict.mexico.1
ggplot(predict.mexico.1, aes(fill=response.level, as.factor(x), predicted)) + 
  ggtitle("Mexico 2012 Panel data") +
  geom_bar(position="stack", stat="identity") + 
  scale_fill_grey(start=0.8, end=0.2, name = " ''I want you to tell me if you agree \n completely, agree somewhat, disagree \n somewhat, or disagree completely: \n My vote is always secret, unless \n I tell someone'' ", labels =  c("Agree completely", "Agree somewhat", "Disagree somewhat", "Disagree completely")) + 
  ylab("Predicted proportion") + xlab("Private transfer offered") + theme_bw() +
  coord_cartesian(expand = TRUE, ylim=c(0, 1))  +
  scale_x_discrete(breaks = c(0,1), labels = c("No", "Yes")) +
  theme(plot.title = element_text(size=12, hjust=0.5), 
        strip.text.x = element_text(size = 9),
        axis.text.x = element_text(size=9), 
        legend.position="left") 
ggsave("mexico2012figure1.png", width=7, height=4, scale=1)



##########################################
##########################################
###   Afrobarometer wave 5 data      ####
##########################################
##########################################
afro5.dat <- read.csv("merged_r5_data_0.csv", header=TRUE, stringsAsFactors=FALSE, row.names=NULL)

afro5.dat$COUNTRY_ALPHA [afro5.dat$COUNTRY_ALPHA==34] <- 10  
afro5.dat$COUNTRY_ALPHA [afro5.dat$COUNTRY_ALPHA==35] <- 29  

afro5.dat$monitored<-NA
afro5.dat$monitored[afro5.dat$Q55==0]<-0
afro5.dat$monitored[afro5.dat$Q55==1]<-1
afro5.dat$monitored[afro5.dat$Q55==2]<-2
afro5.dat$monitored[afro5.dat$Q55==3]<-3

afro5.dat$offer<-NA
afro5.dat$offer[afro5.dat$Q61F==0]<-0
afro5.dat$offer[afro5.dat$Q61F==1]<-1
afro5.dat$offer[afro5.dat$Q61F==2]<-2
afro5.dat$offer[afro5.dat$Q61F==3]<-3

afro5.dat$binary_offer<-NA
afro5.dat$binary_offer[afro5.dat$offer==0]<-0
afro5.dat$binary_offer[afro5.dat$offer>0]<-1

afro5.dat$poverty_food<-NA
afro5.dat$poverty_food[afro5.dat$Q8A==0]<-0
afro5.dat$poverty_food[afro5.dat$Q8A==1]<-1
afro5.dat$poverty_food[afro5.dat$Q8A==2]<-2
afro5.dat$poverty_food[afro5.dat$Q8A==3]<-3
afro5.dat$poverty_food[afro5.dat$Q8A==4]<-4

afro5.dat$poverty_water<-NA
afro5.dat$poverty_water[afro5.dat$Q8B==0]<-0
afro5.dat$poverty_water[afro5.dat$Q8B==1]<-1
afro5.dat$poverty_water[afro5.dat$Q8B==2]<-2
afro5.dat$poverty_water[afro5.dat$Q8B==3]<-3
afro5.dat$poverty_water[afro5.dat$Q8B==4]<-4

afro5.dat$poverty_medicine<-NA
afro5.dat$poverty_medicine[afro5.dat$Q8C==0]<-0
afro5.dat$poverty_medicine[afro5.dat$Q8C==1]<-1
afro5.dat$poverty_medicine[afro5.dat$Q8C==2]<-2
afro5.dat$poverty_medicine[afro5.dat$Q8C==3]<-3
afro5.dat$poverty_medicine[afro5.dat$Q8C==4]<-4

afro5.dat$poverty_fuel<-NA
afro5.dat$poverty_fuel[afro5.dat$Q8D==0]<-0
afro5.dat$poverty_fuel[afro5.dat$Q8D==1]<-1
afro5.dat$poverty_fuel[afro5.dat$Q8D==2]<-2
afro5.dat$poverty_fuel[afro5.dat$Q8D==3]<-3
afro5.dat$poverty_fuel[afro5.dat$Q8D==4]<-4

afro5.dat$poverty_cash<-NA
afro5.dat$poverty_cash[afro5.dat$Q8E==0]<-0
afro5.dat$poverty_cash[afro5.dat$Q8E==1]<-1
afro5.dat$poverty_cash[afro5.dat$Q8E==2]<-2
afro5.dat$poverty_cash[afro5.dat$Q8E==3]<-3
afro5.dat$poverty_cash[afro5.dat$Q8E==4]<-4

afro5.dat$poverty<-(afro5.dat$poverty_food+afro5.dat$poverty_water+afro5.dat$poverty_medicine+afro5.dat$poverty_fuel+afro5.dat$poverty_cash)/20

afro5.dat$age<-afro5.dat$Q1
afro5.dat$age[afro5.dat$Q1==998]<-NA
afro5.dat$age[afro5.dat$Q1==999]<-NA
afro5.dat$age[afro5.dat$Q1==-1]<-NA

afro5.dat$female<-afro5.dat$Q101-1

afro5.dat$education<-NA
afro5.dat$education[afro5.dat$Q97==0]<-0
afro5.dat$education[afro5.dat$Q97==1]<-1
afro5.dat$education[afro5.dat$Q97==2]<-2
afro5.dat$education[afro5.dat$Q97==3]<-3
afro5.dat$education[afro5.dat$Q97==4]<-4
afro5.dat$education[afro5.dat$Q97==5]<-5
afro5.dat$education[afro5.dat$Q97==6]<-6
afro5.dat$education[afro5.dat$Q97==7]<-7
afro5.dat$education[afro5.dat$Q97==8]<-8
afro5.dat$education[afro5.dat$Q97==9]<-9

afro5.dat$contact_party<-NA
afro5.dat$contact_party[afro5.dat$Q30D==0]<-0
afro5.dat$contact_party[afro5.dat$Q30D==1]<-1
afro5.dat$contact_party[afro5.dat$Q30D==2]<-2
afro5.dat$contact_party[afro5.dat$Q30D==3]<-3

afro5.dat$partisanship<-NA
afro5.dat$partisanship[afro5.dat$Q89A==0]<-0
afro5.dat$partisanship[afro5.dat$Q89A==1]<-1

afro5.dat$urban<-NA
afro5.dat$urban[afro5.dat$URBRUR==2]<-0
afro5.dat$urban[afro5.dat$URBRUR==1]<-1

afro5.dat$employment<-NA
afro5.dat$employment[afro5.dat$Q96==0]<-0
afro5.dat$employment[afro5.dat$Q96==1]<-0
afro5.dat$employment[afro5.dat$Q96==2]<-1
afro5.dat$employment[afro5.dat$Q96==3]<-1

afro5.dat$participate_rally<-NA
afro5.dat$participate_rally[afro5.dat$Q29A==0]<-0
afro5.dat$participate_rally[afro5.dat$Q29A==1]<-1

afro5.dat$eval_democracy <-NA
afro5.dat$eval_democracy[afro5.dat$Q42==1] <-0
afro5.dat$eval_democracy[afro5.dat$Q42==2] <-1
afro5.dat$eval_democracy[afro5.dat$Q42==3] <-2
afro5.dat$eval_democracy[afro5.dat$Q42==4] <-3

afro5.dat$turnout<-NA
afro5.dat$turnout[afro5.dat$Q27==0| afro5.dat$Q27==2 | afro5.dat$Q27==3 | afro5.dat$Q27==4 | afro5.dat$Q27==5 | afro5.dat$Q27==6 | afro5.dat$Q27==7 | afro5.dat$Q27==8] <-0
afro5.dat$turnout[afro5.dat$Q27==1] <-1


######################
###  Appendix B.2 ###
#######################
temp.afro.data <- afro5.dat[c("monitored", "offer","binary_offer",  "poverty","age", "female", "education", "employment", "urban", "eval_democracy", "turnout", "participate_rally", "contact_party", "partisanship")]

sink("table_Appendix_B_2.tex")
stargazer(temp.afro.data, 
          type = "latex",
          covariate.labels = c("Belief of vote choice being monitored", "Frequency of private transfers being offered","Binary indicator of a private transfer being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship"), 
          summary.stat = c("n", "mean", "sd", "min", "max"),
          title="Descriptive statistics of Afro Barometer Wave 5 data", 
          digits=3 
)
sink()

######################
###  Appendix B.3 ###
#######################
afro5.reg.1 <- polr(as.factor(monitored) ~  offer
                    + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.2 <- polr(as.factor(monitored) ~ offer
                    + poverty + age + female + education + employment + urban
                    + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.3 <- polr(as.factor(monitored) ~ offer 
                    + poverty + age + female + education + employment + urban
                    + eval_democracy + turnout + participate_rally + contact_party + partisanship
                    + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

sink("table_Appendix_B_3.tex")
stargazer(afro5.reg.1, afro5.reg.2, afro5.reg.3, 
          type = "latex",
          title="Regression results using Afro Barometer Wave 5 data", 
          omit = "COUNTRY_ALPHA",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("Country dummies", "Included", "Included", "Included", "Included", "Included", "Included")),
          covariate.labels = c("Frequency of private transfers being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship", "Binary indicator of a private transfer being offered"), 
          #notes        = "All models use ordered logistic regression; cutoff reports are suppressed * p$<$0.05, ** p$<$0.01", 
          notes.append = FALSE, 
          notes.align = "l")

predict.afro.1 <- ggpredict(afro5.reg.3, terms=c("offer"))
predict.afro.1
sink()

###############################
###  Figure 1. Lower panel ###
###############################
ggplot(predict.afro.1, aes(fill=response.level, as.factor(x), predicted)) + 
  ggtitle("Afrobarometer Wave 5 data") +
  geom_bar(position="stack", stat="identity") + 
  scale_fill_grey(start=0.8, end=0.2, name = " ''How likely do you think it is \n that powerful people can \n find out how you voted, even \n though  there is supposed to be \n a secret ballot in this country?'' ", labels =  c("Not at all likely ", "Not very likely", "Somewhat likely", "Very likely")) + 
  ylab("Predicted proportion") + xlab("Frequency of private transfers being offered") + theme_bw() +
  coord_cartesian(expand = TRUE, ylim=c(0, 1))  +
  scale_x_discrete(breaks = c(0,1, 2, 3), labels = c("Never","Once or twice", "A few times", "Often")) + 
  theme(plot.title = element_text(size=12, hjust=0.5), 
        strip.text.x = element_text(size = 9),
        axis.text.x = element_text(size=9)) 

ggsave("afrofigure1.png", width=7, height=4, scale=1)


#######################
### Appendix B.4  #####
#######################
## for the second table (on corruption)
afro5.dat$belief_corrupt_president <-afro5.dat$Q60A
afro5.dat$belief_corrupt_president [ afro5.dat$Q60A < 0] <- NA
afro5.dat$belief_corrupt_president [ afro5.dat$Q60A > 3] <- NA

afro5.dat$belief_corrupt_parliament <-afro5.dat$Q60B
afro5.dat$belief_corrupt_parliament [ afro5.dat$Q60B < 0] <- NA
afro5.dat$belief_corrupt_parliament [ afro5.dat$Q60B > 3] <- NA

afro5.dat$belief_corrupt_gov_official <-afro5.dat$Q60C
afro5.dat$belief_corrupt_gov_official [ afro5.dat$Q60C < 0] <- NA
afro5.dat$belief_corrupt_gov_official [ afro5.dat$Q60C > 3] <- NA

afro5.dat$belief_corrupt_police <-afro5.dat$Q60E
afro5.dat$belief_corrupt_police [ afro5.dat$Q60E < 0] <- NA
afro5.dat$belief_corrupt_police [ afro5.dat$Q60E > 3] <- NA

afro5.dat$belief_corrupt_tax_official <-afro5.dat$Q60F
afro5.dat$belief_corrupt_tax_official [ afro5.dat$Q60F < 0] <- NA
afro5.dat$belief_corrupt_tax_official [ afro5.dat$Q60F > 3] <- NA

afro5.dat$belief_corrupt_judge <-afro5.dat$Q60G
afro5.dat$belief_corrupt_judge [ afro5.dat$Q60G < 0] <- NA
afro5.dat$belief_corrupt_judge [ afro5.dat$Q60G > 3] <- NA

afro5.reg.extra.1 <- polr(as.factor(belief_corrupt_president) ~ offer 
                    + poverty + age + female + education + employment + urban
                    + eval_democracy + turnout + participate_rally + contact_party + partisanship
                    + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.2 <- polr(as.factor(belief_corrupt_parliament) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.3 <- polr(as.factor(belief_corrupt_gov_official) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.4 <- polr(as.factor(belief_corrupt_police) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.5 <- polr(as.factor(belief_corrupt_tax_official) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.6 <- polr(as.factor(belief_corrupt_judge) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

#Six models are divided for running stargazer. 
sink("table_Appendix_B_4_1.tex")
stargazer(afro5.reg.extra.1, afro5.reg.extra.2, afro5.reg.extra.3, 
          type = "latex",
          title="Regression results using Afro Barometer Wave 5 data", 
          omit = "COUNTRY_ALPHA",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("Country dummies", "Included", "Included", "Included", "Included", "Included", "Included")),
          covariate.labels = c("Frequency of private transfers being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship", "Binary indicator of a private transfer being offered"), 
          notes.append = FALSE, 
          notes.align = "l")

stargazer(afro5.reg.extra.4, afro5.reg.extra.5, afro5.reg.extra.6, 
          type = "latex",
          title="Regression results using Afro Barometer Wave 5 data", 
          omit = "COUNTRY_ALPHA",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("Country dummies", "Included", "Included", "Included", "Included", "Included", "Included")),
          covariate.labels = c("Frequency of private transfers being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship", "Binary indicator of a private transfer being offered"), 
          notes.append = FALSE, 
          notes.align = "l")
sink()


## for the second table (on trust)
afro5.dat$trust_president <-afro5.dat$Q59A
afro5.dat$trust_president [ afro5.dat$Q59A < 0] <- NA
afro5.dat$trust_president [ afro5.dat$Q59A > 3] <- NA

afro5.dat$trust_parliament <-afro5.dat$Q59B
afro5.dat$trust_parliament [ afro5.dat$Q59B < 0] <- NA
afro5.dat$trust_parliament [ afro5.dat$Q59B > 3] <- NA

afro5.dat$trust_electoral_commission <-afro5.dat$Q59C
afro5.dat$trust_electoral_commission [ afro5.dat$Q59C < 0] <- NA
afro5.dat$trust_electoral_commission [ afro5.dat$Q59C > 3] <- NA

afro5.dat$trust_tax_department <-afro5.dat$Q59D
afro5.dat$trust_tax_department [ afro5.dat$Q59D < 0] <- NA
afro5.dat$trust_tax_department [ afro5.dat$Q59D > 3] <- NA

afro5.dat$trust_police <-afro5.dat$Q59H
afro5.dat$trust_police [ afro5.dat$Q59H < 0] <- NA
afro5.dat$trust_police [ afro5.dat$Q59H > 3] <- NA

afro5.dat$trust_court <-afro5.dat$Q59J
afro5.dat$trust_court [ afro5.dat$Q59J < 0] <- NA
afro5.dat$trust_court [ afro5.dat$Q59J > 3] <- NA

afro5.reg.extra.7 <- polr(as.factor(trust_president) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.8 <- polr(as.factor(trust_parliament) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.9 <- polr(as.factor(trust_electoral_commission) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.10 <- polr(as.factor(trust_tax_department) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.11 <- polr(as.factor(trust_police) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

afro5.reg.extra.12 <- polr(as.factor(trust_court) ~ offer 
                          + poverty + age + female + education + employment + urban
                          + eval_democracy + turnout + participate_rally + contact_party + partisanship
                          + as.factor(COUNTRY_ALPHA), data = afro5.dat, Hess = "TRUE")

sink("table_Appendix_B_4_2.tex")
stargazer(afro5.reg.extra.7, afro5.reg.extra.8, afro5.reg.extra.9, 
          type = "latex",
          title="Regression results using Afro Barometer Wave 5 data", 
          omit = "COUNTRY_ALPHA",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("Country dummies", "Included", "Included", "Included", "Included", "Included", "Included")),
          covariate.labels = c("Frequency of private transfers being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship", "Binary indicator of a private transfer being offered"), 
          notes.append = FALSE, 
          notes.align = "l")

stargazer(afro5.reg.extra.10, afro5.reg.extra.11, afro5.reg.extra.12, 
          type = "latex",
          title="Regression results using Afro Barometer Wave 5 data", 
          omit = "COUNTRY_ALPHA",
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("Country dummies", "Included", "Included", "Included", "Included", "Included", "Included")),
          covariate.labels = c("Frequency of private transfers being offered", "Poverty", "Age",  "Female", "Education", "Employment", "Urban", "Evaluation of the country's democracy", "Turnout", "Participation in electoral campaign", "Contact with party", "Partisanship", "Binary indicator of a private transfer being offered"), 
          notes.append = FALSE, 
          notes.align = "l")
sink()

################################################
################################################
###### Mexico survey experiment  ################
################################################
################################################
data<-read.csv("Mexico-Survey-2019_June.csv", header=TRUE, stringsAsFactors=FALSE, row.names=NULL)

data$age<-NA
for (i in 1:length(data$Q1)){
  if (data$Q1[i]==1){
    if (is.na(data$Q1.1[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1[i]}
  }
  if (data$Q1[i]==2){
    if (is.na(data$Q1.1.1[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1.1[i]}
  }
  if (data$Q1[i]==3){
    if (is.na(data$Q1.1.2[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1.2[i]}
  }
  if (data$Q1[i]==4){
    if (is.na(data$Q1.1.3[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1.3[i]}
  }
  if (data$Q1[i]==5){
    if (is.na(data$Q1.1.4[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1.4[i]}
  }
  if (data$Q1[i]==6){
    if (is.na(data$Q1.1.5[i])){data$age[i]<-data$age[i]}
    else {data$age[i]<-data$Q1.1.5[i]}
  }
}

data$state<-data$Q2

data$male<-data$Q4
data$male[data$Q4==99]<-NA

data$education<-data$Q5
data$education[data$Q5==99]<-NA
data$education <- data$education -1

data$full_time_job<-NA
data$full_time_job[data$Q6==0 | data$Q6==1 |data$Q6==3 ] <-0
data$full_time_job[data$Q6==2 ] <-1

data$part_time_job<-NA
data$part_time_job[data$Q6==0 | data$Q6==2 |data$Q6==3 ] <-0
data$part_time_job[data$Q6==1 ] <-1

data$any_job <-NA
data$any_job[data$Q6==0 | data$Q6==3 ]<-0
data$any_job[data$Q6==1 | data$Q6==2 ]<-1

data$positive_reciprocity<-data$Q11
data$positive_reciprocity[data$Q11==99] <-NA

data$negative_reciprocity<-data$Q12
data$negative_reciprocity[data$Q12==99] <-NA
data$negative_reciprocity <- 4- data$negative_reciprocity

data$income<-data$Q7
data$income[data$Q7==99] <-NA

data$any_religion<-NA
data$any_religion[data$Q8==0]<-0
data$any_religion[data$Q8>0 & data$Q8<99]<-1

data$catholic<-NA
data$catholic[data$Q8==0 | data$Q8==2 | data$Q8==3 | data$Q8==4 | data$Q8==5 | data$Q8==6 | data$Q8==7]<-0
data$catholic[data$Q8==1] <-1

data$trust_electoralcommission<-data$Q9.2_2
data$trust_electoralcommission[data$Q9.2_2==99] <-NA

data$turnout_2018_election<-data$Q14
data$turnout_2018_election[data$Q14==1] <-0
data$turnout_2018_election[data$Q14==2] <-1
data$turnout_2018_election[data$Q14==99] <-NA

data$belief_vote_value <-data$Q18
data$belief_vote_value[data$Q18==99] <-NA
data$belief_vote_value <- 4-data$belief_vote_value 

data$offered_2018_election <-data$Q21.2
data$offered_2018_election[data$Q21.2==88 | data$Q21.2==99] <-NA

data$belief_monitored_control <- data$Q24.A
data$belief_monitored_control[data$Q24.A ==99] <-NA

data$belief_monitored_manipulation<- data$Q24.B
data$belief_monitored_manipulation[data$Q24.B ==99] <-NA

data$exp_condition <-NA
data$exp_condition[data$Q24.A >=0] <-0
data$exp_condition[data$Q24.B >=0] <-1

data$belief_monitored <- NA
for (i in 1:nrow(data)){
  if (data$exp_condition[i]==0){data$belief_monitored[i] <- data$belief_monitored_control[i]}
  else {data$belief_monitored[i] <- data$belief_monitored_manipulation[i]}
}

data$correct_manipulation_check <-0 ## refuse to answer is coded as "incorrect"
for (i in 1:nrow(data)){
  if (data$exp_condition[i]==0){
    if (data$Q25[i]==0){
      data$correct_manipulation_check[i]<-1
    }
  }
  else {
    if (data$Q25[i]==1){
      data$correct_manipulation_check[i]<-1
    }
  }
}

#######################
### Appendix C.3  #####
#######################
temp.mexico.stagazer <- data[c("exp_condition", "belief_monitored", "age","male", "education", "income", "any_religion", "any_job","positive_reciprocity", "negative_reciprocity", "belief_vote_value", "trust_electoralcommission", "turnout_2018_election", "offered_2018_election")]

sink("table_Appendix_C_3.tex")
stargazer(temp.mexico.stagazer, 
          type = "latex",
          covariate.labels = c("Manipulation", "Belief of vote choice being monitored", "Age", "Male", "Education", "Income", "Any religion", "Employment status", 
                               "Positive reciprocity", "Negative reciprocity", "Confidence in the impact of my vote",
                               "Trust in electoral commission", "Turnout in the 2018 election", "Offered in the 2018 election"),
          title="Descriptive statistics of Mexico survey data", 
          summary.stat = c("n", "mean", "sd", "min", "max"),
          digits=3 
)
sink()

#######################
### Appendix C.4  #####
#######################
balance.table<-as.data.frame(matrix(ncol = 3, nrow = 10))

ttest.var<- c("age", "male", "education", "income", "any_religion", "any_job", "positive_reciprocity", "negative_reciprocity", "belief_vote_value", "trust_electoralcommission", "turnout_2018_election", "offered_2018_election")

for (i in 1:length(ttest.var)){
  temp.var <- paste("data","$", ttest.var[i], sep="")
  ttest.outcome <- t.test(eval(parse(text=temp.var)) ~ data$exp_condition)
  balance.table[i, 1] <- round(ttest.outcome$estimate[[1]], 3)
  balance.table[i, 2] <- round(ttest.outcome$estimate[[2]], 3)
  
  if (ttest.outcome$p.value<=0.01){
    balance.table[i, 3] <- paste(round(ttest.outcome$estimate[[1]] - ttest.outcome$estimate[[2]], 3), "**", sep="")
  }
  if (ttest.outcome$p.value>0.01 & ttest.outcome$p.value<=0.05){
    balance.table[i, 3] <- paste(round(ttest.outcome$estimate[[1]] - ttest.outcome$estimate[[2]], 3), "*", sep="")
  }
  if (ttest.outcome$p.value>0.05){
    balance.table[i, 3] <- round(ttest.outcome$estimate[[1]] - ttest.outcome$estimate[[2]], 3)
  }
}

ttest.var.label <- c("Age", "Male", "Education", "Income", "Any religion", "Employment status","Positive reciprocity", "Negative reciprocity", "Confidence in the impact of my vote", "Trust in electoral commission", "Turnout in the 2018 election", "Offered in the 2018 election")
balance.output <-cbind(ttest.var.label, balance.table)

balance<- kable(balance.output, "latex",
                caption = "Balance test (Mexico Survey Experiment)", 
                col.names = c("Variable", "Control Mean", "Manipulation Mean", "Difference in Mean"),
                digits = 3 ) %>% 
  kable_styling(bootstrap_options = "striped", full_width = T) %>% 
  footnote(general = "* 0.05, ** 0.01", general_title = "T-test used") 
sink("table_Appendix_C_4.tex")
print(balance)
sink()

###################
### Figure 2  #####
####################
histo_data<-as.data.frame(matrix(ncol =3, nrow = 8))
colnames(histo_data)<-c("belief_monitored", "exp_condition", "value")
histo_data$exp_condition<-rep(c(0, 1), 4)
histo_data$belief_monitored<-c(0, 0, 1, 1, 2, 2, 3, 3)

total_baseline <- nrow(subset(data[c("belief_monitored", "exp_condition")], exp_condition==0 & belief_monitored >=0))
total_manipulation <- nrow(subset(data[c("belief_monitored", "exp_condition")], exp_condition==1 & belief_monitored >=0))

histo_data$value[1]<-nrow(subset(data, exp_condition==0 & belief_monitored==0)) /total_baseline
histo_data$value[2]<-nrow(subset(data, exp_condition==1 & belief_monitored==0)) /total_manipulation 
histo_data$value[3]<-nrow(subset(data, exp_condition==0 & belief_monitored==1)) /total_baseline
histo_data$value[4]<-nrow(subset(data, exp_condition==1 & belief_monitored==1)) /total_manipulation 
histo_data$value[5]<-nrow(subset(data, exp_condition==0 & belief_monitored==2)) /total_baseline
histo_data$value[6]<-nrow(subset(data, exp_condition==1 & belief_monitored==2)) /total_manipulation 
histo_data$value[7]<-nrow(subset(data, exp_condition==0 & belief_monitored==3)) /total_baseline
histo_data$value[8]<-nrow(subset(data, exp_condition==1 & belief_monitored==3)) /total_manipulation 

histo<-ggplot(histo_data, aes(fill=factor(exp_condition), y=value, x=factor(belief_monitored))) + 
  geom_bar(position=position_dodge(), stat="identity") + theme_bw() +
  scale_fill_manual(values=c("#999999", "#666666"), name = "Experimental Condition", labels = c("Baseline (Not offered, 516 obs.)", "Manipulation (Offered, 514 obs.)")) +
  theme(legend.position = c(0.7, 0.85), legend.title = element_text(size = 12), legend.text = element_text(size = 10), 
        axis.text.x = element_text(size=12), axis.text.y = element_text(size=10), axis.title.y = element_text(size = 12)) +
  scale_x_discrete(name = "", labels=c("0" = "Not likely at all", "1" = "A little likely", "2" = "Somewhat likely", "3" = "Very likely")) +
  scale_y_continuous(name="Proportion in each condition", limits=c(0, 0.5)) 

temp.data <- data %>% select(exp_condition, belief_monitored) 
my_sum <- temp.data %>%
  group_by(exp_condition) %>%
  summarise( 
    n=n(),
    mean=mean(belief_monitored, na.rm = TRUE),
    sd=sd(belief_monitored, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
my_sum$color <- my_sum$exp_condition

mean.bar <- ggplot(my_sum, aes(x=factor(exp_condition), y=mean, fill = factor(color))) +
  geom_bar(stat="identity") + theme_bw() +
  geom_errorbar( aes(x=factor(exp_condition), ymin=mean-ic, ymax=mean+ic), width=0.2, size=0.4) +
  scale_fill_manual(values=c("#999999", "#666666"), name = "Experimental Condition", labels = c("Baseline (Not offered, 516 obs.)", "Manipulation (Offered, 514 obs.)")) +
  theme(legend.position = "none", axis.text.x = element_text(size=12), axis.text.y = element_text(size=10), axis.title.y = element_text(size = 12)) +
  scale_x_discrete(name="", limits=factor(c(0, 1)), labels=c("0" = "Baseline", "1" = "Manipulation")) + 
  scale_y_continuous(name="Mean score", limits=c(0, 3), labels=c("0" = "Not likely\n at all", "1" = "A little\n likely", "2" = "Somewhat\n likely", "3" = "Very\n likely")) 

t.test(data$belief_monitored ~ data$exp_condition)

figure <- ggarrange(mean.bar , histo, widths = c(1,2))
annotate_figure(figure, 
                bottom = text_grob(" ``If the broker wanted to find out who you vote for, how likely is it that he can actually find out?'' ", 
                                   color = "black", size =12))
ggsave("histo_mean_bar.png", width=8, height=5, scale=1.1)

###################
### Table 1   #####
####################
regression.1 <- polr(as.factor(belief_monitored) ~ exp_condition + as.factor(state)
                     , data, Hess = "TRUE")

regression.2 <- polr(as.factor(belief_monitored) ~ exp_condition 
                     + age + male + education + income + any_religion + any_job
                     + as.factor(state)
                     , data, Hess = "TRUE")

regression.3 <- polr(as.factor(belief_monitored) ~ exp_condition 
                     + age + male + education + income + any_religion + any_job 
                     + positive_reciprocity + negative_reciprocity
                     + as.factor(state)
                     , data, Hess = "TRUE")

regression.4 <- polr(as.factor(belief_monitored) ~ exp_condition 
                     + age + male + education + income + any_religion + any_job 
                     + positive_reciprocity + negative_reciprocity
                     + belief_vote_value + trust_electoralcommission + turnout_2018_election + offered_2018_election
                     + as.factor(state)
                     , data, Hess = "TRUE")

sink("table_1.tex")
stargazer(regression.1, regression.2, regression.3, regression.4, 
          type = "latex",
          title = "Regression results, using survey experiment data from Mexico", 
          covariate.labels = c("Manipulation vignette",
                               "Age", "Male", "Education", "Income", "Any religion", "Employment status", 
                               "Positive reciprocity", "Negative reciprocity",
                               "Confidence in the impact of my vote", "Trust in electoral commission",
                               "Turnout in the 2018 election", "Offered in the 2018 election"),
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          add.lines = list(c("State dummies", "Included", "Included", "Included", "Included")), 
          omit = "state",
          column.labels = c("Baseline", "Demographic", "Social preference", "Poli. belief and exp."), 
          notes        = "All models use ordered logistic regression; coefficients are non-exponentiated; cutoffs are suppressed in the report; * p$<$0.05, ** p$<$0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()


#############
# FIGURE 3
#############

### Cost is common to Figure 3A and 3B
cost <- 20 # Cost of voting

### Figure 3A
datosA <- data.frame(benefitA = seq(from=0, to=40, length.out = 1000)) # Expressive benefit b_A from 0 to 40
datosA$sanctionModerate <- cost - datosA$benefitA # Sanction for weak supporter of A (b_A < c)
datosA$sanctionStrong <- NA # Sanction for strong supporter of A (b_A > c)
datosA$sanction <- ifelse(datosA$benefitA > cost,
                          datosA$sanctionStrong, datosA$sanctionModerate) # Assign sanctions accordingly

# Create shaded areas (areas above sanction line)
shadedareaA <- data.frame(x=c(0, cost, cost, 0), y=c(cost, 0, 2*cost + 2*cost, 2*cost + 2*cost), type = c('a','a','a','a')) 

# Plot
separatingA <- ggplot(datosA, aes(x = benefitA, y = sanction)) +
  geom_polygon(data = shadedareaA, aes(x = x, y = y), fill = "gray70") +
  geom_line(size = 1.25, col = "black") +
  geom_vline(xintercept = cost,  col = "black", lty = "dashed", size = 0.75) +
  scale_x_continuous(expression(paste("\n Expressive utility from voting for candidate ", italic("A"),", ", italic(b[A]))), limits = c(0, 40), breaks=c(0,20), labels = c(0,expression(paste(bold(italic("c"))))), expand = c(0, 0)) + 
  scale_y_continuous(expression(paste("Sanction from reneging, ", italic("s"), "\n")), breaks=c(0,20), labels = c("",expression(paste(italic("c")))), limits = c(0, 80),  expand = c(0, 0)) + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size=11),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        axis.text.y = element_text(size=13, face = "bold"),
        axis.text.x = element_text(size=13, face = "bold"),
        legend.direction = "horizontal", legend.position = "", legend.title = element_blank()) +
  guides(color = guide_legend(title.position="top", title.hjust=.5, keywidth = 2, keyheight = 1, reverse=FALSE)) 
separatingA

ggsave("separatingA.pdf", width = 5, height = 4.5, unit = "in")

### Figure 3B
datosB <- data.frame(benefitB = seq(from=0, to=40, length.out = 1000)) # Expressive benefit b_B from 0 to 40
datosB$sanctionModerate <- cost + datosB$benefitB # Sanction for weak supporter of B (b_B < c)
datosB$sanctionStrong <- 2 * datosB$benefitB # Sanction for strong supporter of B (b_B > c)
datosB$sanction <- ifelse(datosB$benefitB > cost,
                          datosB$sanctionStrong, datosB$sanctionModerate) # Assign sanctions accordingly

# Create shaded areas (areas above sanction line)
shadedareaB <- data.frame(x=c(0, cost, cost, 0, cost, 2*cost, cost), y=c(cost, 2*cost, 2*cost + 2*cost, 2*cost + 2*cost, 2*cost, 2*cost + 2*cost, 2*cost + 2*cost), type = c('a','a','a','a','b','b','b')) 

# Plot
separatingB <- ggplot(datosB, aes(x = benefitB, y = sanction)) +
  geom_polygon(data = shadedareaB, aes(x = x, y = y, group = as.factor(type), fill = type)) +
  geom_line(size = 1.25, col = "black") +
  geom_vline(xintercept = cost,  col = "black", lty = "dashed", size = 0.75) +
  scale_fill_manual(values = c("gray70", "gray90")) + 
  scale_x_continuous(expression(paste("\n Expressive utility from voting for candidate ", italic("B"),", ",  italic(b[B]))), limits = c(0, 40), breaks=c(0,20), labels = c(0,expression(paste(bold(italic("c"))))), expand = c(0, 0)) + 
  scale_y_continuous(expression(paste("Sanction from reneging, ", italic("s"), "\n")), breaks=c(0,20), labels = c("",expression(paste(italic("c")))), limits = c(0, 80),  expand = c(0, 0)) + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(size=11),
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        axis.text.y = element_text(size=13, face = "bold"),
        axis.text.x = element_text(size=13, face = "bold"),
        legend.direction = "horizontal", legend.position = "", legend.title = element_blank()) +
  guides(color = guide_legend(title.position="top", title.hjust=.5, keywidth = 2, keyheight = 1, reverse=FALSE)) 

separatingB

ggsave("separatingB.pdf", width = 5, height = 4.5, unit = "in")



########################
########################
### Lab experiment  ####
#########################
########################
lab.data <-read.csv("lab_raw_data.csv", header=TRUE, stringsAsFactors=FALSE, row.names=NULL)

lab.data$crt.1 <- 0
lab.data$crt.1[lab.data$crt_answer1==5] <-1
lab.data$crt.2 <- 0
lab.data$crt.2[lab.data$crt_answer2==5] <-1
lab.data$crt.3 <- 0
lab.data$crt.3[lab.data$crt_answer3==47] <-1
lab.data$crt_total_correct <- lab.data$crt.1 + lab.data$crt.2 + lab.data$crt.3

lab.data$female <- lab.data$gender

lab.data$christian <- 0
lab.data$christian[lab.data$religion==1 | lab.data$religion==2 |lab.data$religion==3] <-1

lab.data$muslim <- 0
lab.data$muslim[lab.data$religion ==5] <-1

lab.data$total_election_earning <- lab.data$earning_election1+lab.data$earning_election2+lab.data$earning_election3+lab.data$earning_election4+lab.data$earning_election5+lab.data$earning_election6+lab.data$earning_election7+lab.data$earning_election8+lab.data$earning_election9+lab.data$earning_election10+lab.data$earning_election11+lab.data$earning_election12

## strategy classification
lab.data$strategy <- NA
lab.data$strategy[lab.data$belieftask_belief_notransfer==1 | lab.data$belieftask_belief_transfer==0 ] <- 0
lab.data$strategy[lab.data$belieftask_belief_notransfer==0 & lab.data$belieftask_belief_transfer==1 ] <- 1
lab.data$strategy[lab.data$belieftask_belief_notransfer==0 & lab.data$belieftask_belief_transfer==0 ] <- 2
lab.data$strategy[lab.data$belieftask_belief_notransfer==1 & lab.data$belieftask_belief_transfer==1 ] <- 3

lab.data$odd.strategy <- 0
lab.data$odd.strategy[lab.data$strategy == 0] <-1 
lab.data$odd.strategy[is.na(lab.data$strategy)] <- NA

lab.data$separating <- 0
lab.data$separating[lab.data$strategy == 1] <-1
lab.data$separating[is.na(lab.data$strategy)] <-NA

lab.data$pooling_no_monitoring <- 0
lab.data$pooling_no_monitoring[lab.data$strategy == 2] <-1
lab.data$pooling_no_monitoring[is.na(lab.data$strategy)] <-NA

lab.data$pooling_monitoring <-0
lab.data$pooling_monitoring[lab.data$strategy == 3] <-1
lab.data$pooling_monitoring[is.na(lab.data$strategy)] <-NA


###################
### Figure 4  #####
####################
## left panel
lab.data$belieftask_belief_notransfer[lab.data$belieftask_belief_notransfer==-99] <-NA
lab.data$belieftask_belief_transfer[lab.data$belieftask_belief_transfer==-99] <-NA

dt.1 <- lab.data%>%
  dplyr::summarise(
    mean = mean(belieftask_belief_notransfer, na.rm =TRUE),
    lci = t.test(belieftask_belief_notransfer, conf.level = 0.95)$conf.int[1],
    uci = t.test(belieftask_belief_notransfer, conf.level = 0.95)$conf.int[2])
dt.1

dt.2 <- lab.data%>%
  dplyr::summarise(
    mean = mean(belieftask_belief_transfer, na.rm =TRUE),
    lci = t.test(belieftask_belief_transfer, conf.level = 0.95)$conf.int[1],
    uci = t.test(belieftask_belief_transfer, conf.level = 0.95)$conf.int[2])
dt.2

temp.belief.data.1 <- data.frame(
  Transfer = c("Not offered", "Offered"), 
  Proportion = c(dt.1[1,1], dt.2[1,1]), 
  lci = c(dt.1[1,2], dt.2[1,2]), 
  uci = c(dt.1[1,3], dt.2[1,3])
)
temp.belief.data.1

belief.distribution.1<- ggplot(temp.belief.data.1, aes(x=Transfer, y=Proportion, fill=c("Not offered", "Offered"))) +  
  ggtitle("Binary measure (0/1)") + 
  geom_bar(stat = "identity") + theme_bw() + 
  geom_errorbar(aes(x=Transfer, ymin=lci, ymax= uci), width = 0.2, size = 0.5) + 
  scale_y_continuous(name="Proportion", limits=c(0, 1)) +
  scale_fill_manual(values =c("#999999", "#666666")) +
  theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=11), axis.title.y = element_text(size = 12)) +
  guides(fill=FALSE)

# right panel
lab.data$belieftask_confidence_notransfer[lab.data$belieftask_confidence_notransfer==-99] <-NA
lab.data$belieftask_confidence_transfer[lab.data$belieftask_confidence_transfer==-99] <-NA

lab.data$belief_continous_notransfer <- NA
lab.data$belief_continous_transfer <- NA

for (i in 1:nrow(lab.data)){
  if (is.na(lab.data$belieftask_belief_notransfer[i])){
    lab.data$belief_continous_notransfer[i] <-NA
  }
  else {
    if (lab.data$belieftask_belief_notransfer[i]==0){
      lab.data$belief_continous_notransfer[i] <- (5*(10-lab.data$belieftask_confidence_notransfer[i]))/100
    }
    if  (lab.data$belieftask_belief_notransfer[i]==1){
      lab.data$belief_continous_notransfer[i] <- (100 - 5*(10-lab.data$belieftask_confidence_notransfer[i]))/100
    }
  }
}

for (i in 1:nrow(lab.data)){
  if (is.na(lab.data$belieftask_belief_transfer[i])){
    lab.data$belief_continous_transfer[i] <-NA
  }
  else {
    if (lab.data$belieftask_belief_transfer[i]==0){
      lab.data$belief_continous_transfer[i] <- (5*(10-lab.data$belieftask_confidence_transfer[i]))/100
    }
    if  (lab.data$belieftask_belief_transfer[i]==1){
      lab.data$belief_continous_transfer[i] <- (100 - 5*(10-lab.data$belieftask_confidence_transfer[i]))/100
    }
  }
}

dt.3 <- lab.data%>%
  dplyr::summarise(
    mean = mean(belief_continous_notransfer, na.rm =TRUE),
    lci = t.test(belief_continous_notransfer, conf.level = 0.95)$conf.int[1],
    uci = t.test(belief_continous_notransfer, conf.level = 0.95)$conf.int[2])
dt.3

dt.4 <- lab.data%>%
  dplyr::summarise(
    mean = mean(belief_continous_transfer, na.rm =TRUE),
    lci = t.test(belief_continous_transfer, conf.level = 0.95)$conf.int[1],
    uci = t.test(belief_continous_transfer, conf.level = 0.95)$conf.int[2])
dt.4

temp.belief.data.2 <- data.frame(
  Transfer = c("Not offered", "Offered"), 
  Score = c(dt.3[1,1], dt.4[1,1]), 
  lci = c(dt.3[1,2], dt.4[1,2]), 
  uci = c(dt.3[1,3], dt.4[1,3])
)
temp.belief.data.2

belief.distribution.2<- ggplot(temp.belief.data.2, aes(x=Transfer, y=Score, fill=c("Not offered", "Offered"))) +  
  ggtitle("Continuous measure(0~1)") + 
  geom_bar(stat = "identity") + theme_bw() + 
  geom_errorbar(aes(x=Transfer, ymin=lci, ymax= uci), width = 0.2, size = 0.5) + 
  scale_y_continuous(name="Average of belief scores", limits=c(0, 1)) +
  scale_fill_manual(values =c("#999999", "#666666")) +
  theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=11), axis.title.y = element_text(size = 12)) +
  guides(fill=FALSE)

belief_figure <- ggarrange(belief.distribution.1, belief.distribution.2)
belief_figure
ggsave("belieffigure.png", width=7, height=4, scale=1.2)



########################
### Appendix E.2  #####
########################
lab.data$alt_belief_continous_notransfer <- NA
lab.data$alt_belief_continous_transfer <- NA

for (i in 1:nrow(lab.data)){
  if (is.na(lab.data$belieftask_belief_notransfer[i])){
    lab.data$alt_belief_continous_notransfer[i] <-NA
  }
  else {
    if (lab.data$belieftask_belief_notransfer[i]==0){
      lab.data$alt_belief_continous_notransfer[i] <- (4*(10-lab.data$belieftask_confidence_notransfer[i]))/100
    }
    if  (lab.data$belieftask_belief_notransfer[i]==1){
      lab.data$alt_belief_continous_notransfer[i] <- (100 - 4*(10-lab.data$belieftask_confidence_notransfer[i]))/100
    }
  }
}

for (i in 1:nrow(lab.data)){
  if (is.na(lab.data$belieftask_belief_transfer[i])){
    lab.data$alt_belief_continous_transfer[i] <-NA
  }
  else {
    if (lab.data$belieftask_belief_transfer[i]==0){
      lab.data$alt_belief_continous_transfer[i] <- (4*(10-lab.data$belieftask_confidence_transfer[i]))/100
    }
    if  (lab.data$belieftask_belief_transfer[i]==1){
      lab.data$alt_belief_continous_transfer[i] <- (100 - 4*(10-lab.data$belieftask_confidence_transfer[i]))/100
    }
  }
}

temp.belief.data.2.alt <- data.frame(
  Transfer = c("Not offered", "Offered"), 
  Score= c(mean(lab.data$alt_belief_continous_notransfer, na.rm =TRUE), mean(lab.data$alt_belief_continous_transfer, na.rm =TRUE)) 
)

belief.distribution.2.alt<- ggplot(temp.belief.data.2.alt, aes(x=Transfer, y=Score, fill=c("Not offered", "Offered"))) +  
  geom_bar(stat = "identity") + theme_bw() + 
  scale_y_continuous(name="Average of belief scores (alternative continuous measure)", limits=c(0, 1)) +
  scale_fill_manual(values =c("#999999", "#666666")) +
  theme(axis.text.x = element_text(size=12), axis.text.y = element_text(size=11), axis.title.y = element_text(size = 12)) +
  guides(fill=FALSE)
belief.distribution.2.alt
ggsave("belieffigurealt.png", width=3.5, height=4, scale=1.2)


########################
### Appendix E.4  #####
########################
temp.lab.data.stargazer <- lab.data[c("age","female", "christian", "muslim", "beautytask_choice", "crt_total_correct", "altruism", "positive_reciprocity")]

sink("table_Appendix_E_4.tex")
stargazer(temp.lab.data.stargazer, 
          type = "latex",
          covariate.labels = c("Age", "Female", "Christian", "Muslim", "Choice in beauty contest", "Num. of correct answers in the CRT", "Altruism", "Positive reciprocity"),
          title="Descriptive statistics of lab experiment data", 
          summary.stat = c("n", "mean", "sd", "min", "max"),
          digits=3 
)
sink()

###################
### Table 2  #####
####################
lab.data.belief<-reshape(lab.data, direction='long',
                         varying =list(c(which( colnames(lab.data)=="belieftask_belief_notransfer"): which(colnames(lab.data)=="belieftask_belief_transfer")), 
                                       c(which( colnames(lab.data)=="belieftask_confidence_notransfer"), which( colnames(lab.data)=="belieftask_confidence_transfer")),
                                       c(which( colnames(lab.data)=="belief_continous_notransfer"), which( colnames(lab.data)=="belief_continous_transfer")),
                                       c(which( colnames(lab.data)=="alt_belief_continous_notransfer"), which( colnames(lab.data)=="alt_belief_continous_transfer")) ), 
                         v.names= c("binary_belief", "confidence", "continuous_belief", "alt_continuous_belief"), 
                         timevar = "guessing_transfer", 
                         idvar="subject_id")


reg.lab.1 <- glm(as.factor(binary_belief)~ guessing_transfer + role_first_half 
                        , 
                        data = lab.data.belief, 
                        family="binomial")

reg.lab.2 <- glm(as.factor(binary_belief)~ guessing_transfer + role_first_half + age + female + christian + muslim
                        , 
                        data = lab.data.belief, 
                        family="binomial")

reg.lab.3 <- glm(as.factor(binary_belief)~ guessing_transfer + role_first_half + age + female + christian + muslim + beautytask_choice + crt_total_correct + altruism + positive_reciprocity
                        , 
                        data = lab.data.belief, 
                        family="binomial")

reg.lab.4 <- lm(continuous_belief~ guessing_transfer + role_first_half
                       , 
                       data = lab.data.belief)

reg.lab.5 <- lm(continuous_belief~ guessing_transfer + role_first_half + age + female + christian + muslim
                       , 
                       data = lab.data.belief)

reg.lab.6 <- lm(continuous_belief~ guessing_transfer + role_first_half + age + female + christian + muslim + beautytask_choice + crt_total_correct + altruism + positive_reciprocity 
                       , 
                       data = lab.data.belief)

sink("table_2.tex")
stargazer(reg.lab.1, reg.lab.2, reg.lab.3, 
          reg.lab.4, reg.lab.5, reg.lab.6, 
          se = list(sqrt(diag(vcovHC(reg.lab.1 , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.2 , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.3 , type = "HC1"))), 
                    sqrt(diag(vcovHC(reg.lab.4 , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.5 , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.6 , type = "HC1")))),
          type = "latex", 
          title = "Regression results regarding belief Updating",
          covariate.labels = c("Hypothetical scenario of transfer provision", "Role in the first six rounds",
                               "Age", "Female", "Christian", "Muslim", 
                               "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism" ,"Positive reciprocity"),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Binary belief (0/1)", "Continuous belief (0~1)"),
          column.separate = c(3, 3),
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 3 use logit regressions, whereas Models 4 to 6 use OLS regressions; standard errors, clustered at the individual level, in parentheses; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l") 
sink()
ggpredict(reg.lab.3, terms=c("guessing_transfer"))
ggpredict(reg.lab.6, terms=c("guessing_transfer"))


########################
### Appendix E.5  #####
########################
reg.lab.4.alt <- lm(alt_continuous_belief~ guessing_transfer + role_first_half
                           , 
                           data = lab.data.belief)

reg.lab.5.alt <- lm(alt_continuous_belief~ guessing_transfer + role_first_half + age + female + christian + muslim
                           , 
                           data = lab.data.belief)

reg.lab.6.alt <- lm(alt_continuous_belief~ guessing_transfer + role_first_half + age + female + christian + muslim + beautytask_choice + crt_total_correct + altruism + positive_reciprocity 
                           , 
                           data = lab.data.belief)

sink("table_Appendix_E_5.tex")
stargazer(reg.lab.4.alt, reg.lab.5.alt, reg.lab.6.alt, 
          se = list(sqrt(diag(vcovHC(reg.lab.4.alt , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.5.alt , type = "HC1"))), sqrt(diag(vcovHC(reg.lab.6.alt , type = "HC1")))),
          type = "latex", 
          title = "Replication with the alternative measure of a belief along a continuous scale",
          covariate.labels = c("Hypothetical scenario of transfer provision", "Role in the first six rounds",
                               "Age", "Female", "Christian", "Muslim", 
                               "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism" ,"Positive reciprocity"),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Alternative measure of continuous belief (0~1)"),
          column.separate = c(3),
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "All models use OLS regressions; standard errors, clustered at the individual level, in parentheses; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l") 
sink()

#### the followings are to make data with 6 rounds (originally 12 rounds: 2 roles. )
#role_first_half==1: candidate,  role_first_half==2: voter
## candidate decisions of participants of role_first_half==1
lab.data.2<-subset(lab.data, role_first_half==1)
lab.data.3<-reshape(lab.data.2, direction='long',
                    varying = list(c(which( colnames(lab.data)=="transfer1"): which(colnames(lab.data)=="transfer6")),
                                   c(which( colnames(lab.data)=="monitoring_ability1"): which(colnames(lab.data)=="monitoring_ability6")), 
                                   c(which( colnames(lab.data)=="earning_election1"): which(colnames(lab.data)=="earning_election6")), 
                                   c(which( colnames(lab.data)=="election_result1"): which(colnames(lab.data)=="election_result6")), 
                                   c(which( colnames(lab.data)=="partner_vote_choice1"): which(colnames(lab.data)=="partner_vote_choice6")), 
                                   c(which( colnames(lab.data)=="revealed_partner_vote_choice1"): which(colnames(lab.data)=="revealed_partner_vote_choice6")), 
                                   c(which( colnames(lab.data)=="partner_sanctioned1"): which(colnames(lab.data)=="partner_sanctioned6"))),
                    v.names = c("transfer", "monitoring_ability", "earning_election", "election_result", "partner_vote_choice", "revealed_partner_vote_choice", "partner_sanctioned"),
                    timevar = "round",
                    idvar="subject_id")

## candidate decisions of participants of role_first_half==2
lab.data.4<-subset(lab.data, role_first_half==2)
lab.data.5<-reshape(lab.data.4, direction='long',
                    varying = list(c(which( colnames(lab.data)=="transfer7"): which(colnames(lab.data)=="transfer12")),
                                   c(which( colnames(lab.data)=="monitoring_ability7"): which(colnames(lab.data)=="monitoring_ability12")), 
                                   c(which( colnames(lab.data)=="earning_election7"): which(colnames(lab.data)=="earning_election12")), 
                                   c(which( colnames(lab.data)=="election_result7"): which(colnames(lab.data)=="election_result12")), 
                                   c(which( colnames(lab.data)=="partner_vote_choice7"): which(colnames(lab.data)=="partner_vote_choice12")), 
                                   c(which( colnames(lab.data)=="revealed_partner_vote_choice7"): which(colnames(lab.data)=="revealed_partner_vote_choice12")), 
                                   c(which( colnames(lab.data)=="partner_sanctioned7"): which(colnames(lab.data)=="partner_sanctioned12"))),
                    v.names = c("transfer", "monitoring_ability", "earning_election", "election_result", "partner_vote_choice", "revealed_partner_vote_choice", "partner_sanctioned"),
                    timevar = "round",
                    idvar="subject_id")

## vote decisions of participants of role_first_half==1 
lab.data.6<-subset(lab.data, role_first_half==1)
lab.data.7<-reshape(lab.data.6, direction='long',
                    varying = list(c(which( colnames(lab.data)=="vote_choice7"): which(colnames(lab.data)=="vote_choice12")),
                                   c(which( colnames(lab.data)=="sanctioned7"): which(colnames(lab.data)=="sanctioned12")), 
                                   c(which( colnames(lab.data)=="partner_transfer7"): which(colnames(lab.data)=="partner_transfer12")), 
                                   c(which( colnames(lab.data)=="partner_monitoring_ability7"): which(colnames(lab.data)=="partner_monitoring_ability12"))),         
                    v.names = c("vote_choice", "sanctioned", "partner_transfer", "partner_monitoring_ability"),
                    timevar = "round",
                    idvar="subject_id")

lab.data.7_sub<-(lab.data.7[c("subject_id", "round", "vote_choice", "sanctioned", "partner_transfer", "partner_monitoring_ability")])
d3_d7<-merge(lab.data.3, lab.data.7_sub, by = c("subject_id", "round"))
d3_d7[, (which( colnames(d3_d7)=="transfer7"): which(colnames(d3_d7)=="partner_sanctioned12"))]<-NULL  # remove different variables to use rbind later

## vote decisions of participants of role_first_half==2 
lab.data.8<-subset(lab.data, role_first_half==2)
lab.data.9<-reshape(lab.data.8, direction='long',
                    varying = list(c(which( colnames(lab.data)=="vote_choice1"): which(colnames(lab.data)=="vote_choice6")),
                                   c(which( colnames(lab.data)=="sanctioned1"): which(colnames(lab.data)=="sanctioned6")), 
                                   c(which( colnames(lab.data)=="partner_transfer1"): which(colnames(lab.data)=="partner_transfer6")), 
                                   c(which( colnames(lab.data)=="partner_monitoring_ability1"): which(colnames(lab.data)=="partner_monitoring_ability6"))),       
                    v.names = c("vote_choice", "sanctioned", "partner_transfer", "partner_monitoring_ability"),
                    timevar = "round",
                    idvar="subject_id")

lab.data.9_sub<-(lab.data.9[c("subject_id", "round", "vote_choice", "sanctioned", "partner_transfer", "partner_monitoring_ability")])
d5_d9<-merge(lab.data.5, lab.data.9_sub, by = c("subject_id", "round"))
d5_d9[, (which( colnames(d5_d9)=="transfer1"): which(colnames(d5_d9)=="partner_sanctioned6"))]<-NULL # remove different variables to use rbind later

lab.data.6rounds<-rbind(d3_d7, d5_d9)

write.csv(lab.data.6rounds, "lab_6rounds.csv") ## this data will be used to run regressions for Appendix E.6 in Stata.

## lagged variables
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_earning_election = dplyr::lag(earning_election, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_sanctioned = dplyr::lag(sanctioned, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_vote_choice = dplyr::lag(vote_choice, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_transfer = dplyr::lag(transfer, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_partner_transfer = dplyr::lag(partner_transfer, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_partner_sanctioned = dplyr::lag(partner_sanctioned, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_partner_vote_choice = dplyr::lag(partner_vote_choice, n = 1, default = NA))
lab.data.6rounds <- lab.data.6rounds %>% group_by(subject_id) %>% mutate(lagged_monitoring_ability = dplyr::lag(monitoring_ability, n = 1, default = NA))

lab.data.6rounds$binary_diff_belief <- lab.data.6rounds$belieftask_belief_transfer - lab.data.6rounds$belieftask_belief_notransfer
lab.data.6rounds$diff_belief <- lab.data.6rounds$belief_continous_transfer - lab.data.6rounds$belief_continous_notransfer
lab.data.6rounds$alt_diff_belief <- lab.data.6rounds$alt_belief_continous_transfer - lab.data.6rounds$alt_belief_continous_notransfer

lab.data.6rounds$vote_choice <- as.factor(lab.data.6rounds$vote_choice)


#################
### Table 3 ####
#################
reg.lab.7.0 <- clmm(vote_choice ~   
                      partner_transfer + 
                      role_first_half +
                      round +
                      (1|subject_id),
                    data = lab.data.6rounds )

reg.lab.7.1 <- clmm(vote_choice ~   
                      partner_transfer + 
                      role_first_half +
                      age + 
                      female + 
                      christian + 
                      muslim + 
                      round + 
                      (1|subject_id),
                    data = lab.data.6rounds )

reg.lab.7 <- clmm(vote_choice ~   
                    partner_transfer + 
                    role_first_half +
                    age + 
                    female + 
                    christian + 
                    muslim + 
                    beautytask_choice + 
                    crt_total_correct + 
                    altruism + 
                    positive_reciprocity + 
                    round + 
                    (1|subject_id),
                  data = lab.data.6rounds )

k1<-plot_model(reg.lab.7, type = "pred", terms = "partner_transfer")
k1[1]

reg.lab.9.0 <- glmer(transfer ~ 
                       monitoring_ability + 
                       role_first_half +
                       round + 
                       (1|subject_id),
                     family = "binomial", data = lab.data.6rounds, 
                     control = glmerControl(
                       optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.9.1 <- glmer(transfer ~ 
                       monitoring_ability + 
                       role_first_half +
                       age + 
                       female + 
                       christian + 
                       muslim + 
                       round + 
                       (1|subject_id),
                     family = "binomial", data = lab.data.6rounds, 
                     control = glmerControl(
                       optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.9 <- glmer(transfer ~ 
                     monitoring_ability + 
                     role_first_half +
                     age + 
                     female + 
                     christian + 
                     muslim + 
                     beautytask_choice + 
                     crt_total_correct + 
                     altruism + 
                     positive_reciprocity + 
                     round + 
                     (1|subject_id),
                   family = "binomial", data = lab.data.6rounds, 
                   control = glmerControl(
                     optimizer ='optimx', optCtrl=list(method='nlminb')))

k2<-plot_model(reg.lab.9, type = "pred", terms = "monitoring_ability ")
k2[1]

coef.7.0 <- coeftest(reg.lab.7.0)
coef.7.1 <- coeftest(reg.lab.7.1)
coef.7 <- coeftest(reg.lab.7)

sink("table_3.tex")
stargazer(coef.7.0, coef.7.1, coef.7 , 
          reg.lab.9.0, reg.lab.9.1, reg.lab.9, 
          title = "Regression results regarding vote choice and transfer provision",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism", 
                                "Positive reciprocity", "Round"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(3, 3), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 6 use random effect regressions (ologit for Models 1 to 3 and logit for Models 4 to 6); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()

#####################
### Appendix E.6 ####
#####################
# For technical reasons, regressions for Appendix E.6 are performed in Stata. 
# For do and log files, see Appendix E6.do and Appendix E6_log.smcl, respectively. 

########################
### Appendix E.7  #####
########################
lab.data$cont_diff <- NA
for (i in 1:nrow(lab.data)){
  if (!is.na(lab.data$belief_continous_transfer[i])) {
    lab.data$cont_diff[i] <- lab.data$belief_continous_transfer[i] - lab.data$belief_continous_notransfer[i]
  }
}
table(lab.data$cont_diff)

ggplot(lab.data, aes(x=cont_diff)) +  theme_bw() +
  geom_histogram(aes (x=cont_diff), binwidth = 0.05, colour = "black", fill = "grey") +
  geom_density(alpha = 0.5) +
  labs(x = "Belief updating", y = "Frequency") +
  theme(text = element_text(size= 12)) + 
  scale_y_continuous(breaks = seq(0, 16, by = 2)) + 
  scale_x_continuous(breaks = seq(-1, 1, by = 0.2)) +
  ggtitle("Panel A")

ggsave("distribution_diff_belief.png", width=7, height=4, scale=1.1)


lab.data$alt_cont_diff <- NA
for (i in 1:nrow(lab.data)){
  if (!is.na(lab.data$alt_belief_continous_transfer[i])) {
    lab.data$alt_cont_diff[i] <- lab.data$alt_belief_continous_transfer[i] - lab.data$alt_belief_continous_notransfer[i]
  }
}
table(lab.data$alt_cont_diff)

ggplot(lab.data, aes(x=alt_cont_diff)) +  theme_bw() +
  geom_histogram(aes (x=alt_cont_diff), binwidth = 0.05, colour = "black", fill = "grey") +
  geom_density(alpha = 0.5) +
  labs(x = "Belief updating", y = "Frequency") +
  theme(text = element_text(size= 14)) + 
  scale_y_continuous(breaks = seq(0, 16, by = 2)) + 
  scale_x_continuous(breaks = seq(-1, 1, by = 0.2)) + 
  ggtitle("Panel B")

ggsave("distribution_alt_diff_belief.png", width=7, height=4, scale=1.1)


########################
### Appendix E.8  #####
########################
reg.lab.8 <- clmm(vote_choice ~   
                           diff_belief *  partner_transfer + 
                           role_first_half +
                           age + 
                           female + 
                           christian + 
                           muslim + 
                           beautytask_choice + 
                           crt_total_correct + 
                           altruism + 
                           positive_reciprocity + 
                           round + 
                           (1|subject_id),
                         data = lab.data.6rounds )

lab.data.6rounds$transfer <- as.factor(lab.data.6rounds$transfer)

reg.lab.10 <- glmer(transfer ~ 
                            diff_belief * monitoring_ability + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            crt_total_correct + 
                            altruism + 
                            positive_reciprocity + 
                            round + 
                            (1|subject_id),
                          family = "binomial", data = lab.data.6rounds, 
                          control = glmerControl(
                            optimizer ='optimx', optCtrl=list(method='nlminb')))

sink("table_Appendix_E_8.tex")
stargazer(coeftest(reg.lab.8), reg.lab.10, 
          title = "Regression results regarding vote choice and transfer provision according to belief updating (adjusted score)",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Diff. in adjusted beliefs", "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism", 
                                "Positive reciprocity", "Round",
                                "Diff. in adjusted beliefs X Private transfer offered" , 
                                "Diff. in adjusted beliefs X Monitoring capacity"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(1, 1), 
          add.lines = list(c("Observations", "360", "360")), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 and 2 use random effect regressions (ologit for Model 1 and logit for Model 2); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()

###################
### Figure 5  #####
###################
response.level.labs<- c("Voting for Programmatic candidate", "Abstention", "Voting for Clientelistic candidate")
names(response.level.labs) <- c("1", "2", "3")

margin.seg.1 <- data.frame(response.level = 3, x1 = 0, x2= 0, y1=0.06, y2=0.56)
margin.seg.4 <- data.frame(response.level = 3, x1 = 1, x2= 1, y1=0.07, y2=0.83)
margin.seg.7 <- data.frame(response.level = 3, x1 = 0.1, y1=0.33, lab ="A")
margin.seg.8 <- data.frame(response.level = 3, x1 = 0.9, y1=0.46, lab ="B")

p1<- plot_model(reg.lab.8, type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "partner_transfer"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Private transfer",
                axis.title = "Predicted probability", show.legend = TRUE) + theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1,0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") + 
  scale_linetype_manual(values = c("solid", "dotted"), 
                        labels = c("Not offered", "Offered")) +
  geom_segment(data = margin.seg.1, aes(x = x1, y = y1, xend=x2, yend = y2), arrow= arrow(length =unit(0.2, "cm"), ends="both"),
               colour='darkgray',size=1,linetype=1, inherit.aes=FALSE) +
  geom_segment(data = margin.seg.4, aes(x = x1, y = y1, xend=x2, yend = y2), arrow= arrow(length =unit(0.2, "cm"), ends="both"), 
               colour='darkgray',size=1,linetype=1, inherit.aes=FALSE) +
  geom_text(data = margin.seg.7, mapping = aes(x = x1, y = y1, label = lab), size = 4, inherit.aes=FALSE) +
  geom_text(data = margin.seg.8, mapping = aes(x = x1, y = y1, label = lab), size = 4, inherit.aes=FALSE) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = response.level.labs)) + theme(strip.text.x = element_text(size = 10))


p2<- plot_model(reg.lab.10 , type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "monitoring_ability"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Monitoring capacity",
                axis.title = "Predicted probability of offering transfer", show.legend = TRUE) +  theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1, 0, 1)) + xlab("Belief updating") + theme(legend.position="bottom")  +
  scale_linetype_manual(values = c("solid", "dotted"), 
                        labels = c("Weak", "Strong"))

ggarrange(p1, p2, nrow = 2)
ggsave("marginal_effect_belief_updating.png", width=5, height=6, scale=1.6)


#####################
### Appendix E.9  ###
#####################
reg.lab.11 <- clmm(vote_choice ~   
                            binary_diff_belief +  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            crt_total_correct + 
                            altruism + 
                            positive_reciprocity + 
                            round + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.12 <- clmm(vote_choice ~   
                            binary_diff_belief *  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            crt_total_correct + 
                            altruism + 
                            positive_reciprocity + 
                            round + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.13 <- glmer(transfer ~ 
                             binary_diff_belief + monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             crt_total_correct + 
                             altruism + 
                             positive_reciprocity + 
                             round + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.14 <- glmer(transfer ~ 
                             binary_diff_belief * monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             crt_total_correct + 
                             altruism + 
                             positive_reciprocity + 
                             round + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

coef.11 <- coeftest(reg.lab.11)
coef.12 <- coeftest(reg.lab.12)

sink("table_Appendix_E_9.tex")
stargazer(coef.11, coef.12, 
          reg.lab.13, reg.lab.14, 
          title = "Replication of Table 3 using the binary measure of a belief",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Diff. in binary beliefs", "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism", 
                                "Positive reciprocity", "Round",
                                "Diff. in binary beliefs X Transfer offered" , 
                                "Diff. in binary beliefs X Monitoring capacity"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(2, 2), 
          add.lines = list(c("Observations", "360", "360", "360", "360")), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 4 use random effect regressions (ologit for Models 1 and 2 and logit for Models 3 and 4); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()

p3 <- ggpredict(reg.lab.12, terms = c("binary_diff_belief[-1,0,1]", "partner_transfer")) 
ggplot(p3, aes(as.factor(x), predicted, group=group), colour=group) + theme_bw() +
  geom_point(aes(shape=group), size=2, position = position_dodge(width = 0.2)) +  
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high), 
    position = position_dodge(.2), 
    width=0.2, 
    size=0.5) +
  scale_x_discrete(breaks = c(-1, 0, 1), labels = c("-1", "0", "1")) + 
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), labels=c("0", "0.25", "0.5", "0.75",  "1")) + 
  xlab("Belief updating")  + ylab("Predicted probability") + scale_shape_discrete(name="Private transfer", labels=c("Not offered", "Offered")) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = response.level.labs)) + 
  theme(legend.position="bottom", 
        strip.text.x = element_text(size = 10))
ggsave("marginaleffect3.png", width=7, height=4, scale=1.1)

p4 <- ggpredict(reg.lab.14, terms = c("binary_diff_belief[-1,0,1]", "monitoring_ability")) 
ggplot(p4, aes(as.factor(x), predicted, group=group), colour=group) + theme_bw() +
  geom_point(aes(shape=group), size=2, position = position_dodge(width = 0.2)) +  
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high), 
    position = position_dodge(.2), 
    width=0.2, 
    size=0.5) +
  scale_x_discrete(breaks = c(-1, 0, 1), labels = c("-1", "0", "1")) + 
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), labels=c("0", "0.25", "0.5", "0.75",  "1")) + 
  xlab("Belief updating")  + ylab("Predicted probability") + scale_shape_discrete(name="Monitoring capacity", labels=c("No", "Yes")) +
  theme(legend.position="bottom", 
        strip.text.x = element_text(size = 10))
ggsave("marginaleffect4.png", width=7, height=4, scale=1.1)


#####################
### Appendix E.10  ###
#####################
reg.lab.7.alt <- clmm(vote_choice ~   
                               alt_diff_belief +  partner_transfer + 
                               role_first_half +
                               age + 
                               female + 
                               christian + 
                               muslim + 
                               beautytask_choice + 
                               crt_total_correct + 
                               altruism + 
                               positive_reciprocity + 
                               round + 
                               (1|subject_id),
                             data = lab.data.6rounds )

reg.lab.8.alt <- clmm(vote_choice ~   
                               alt_diff_belief *  partner_transfer + 
                               role_first_half +
                               age + 
                               female + 
                               christian + 
                               muslim + 
                               beautytask_choice + 
                               crt_total_correct + 
                               altruism + 
                               positive_reciprocity + 
                               round + 
                               (1|subject_id),
                             data = lab.data.6rounds )

reg.lab.9.alt <- glmer(transfer ~ 
                                alt_diff_belief + monitoring_ability + 
                                role_first_half +
                                age + 
                                female + 
                                christian + 
                                muslim + 
                                beautytask_choice + 
                                crt_total_correct + 
                                altruism + 
                                positive_reciprocity + 
                                round + 
                                (1|subject_id),
                              family = "binomial", data = lab.data.6rounds, 
                              control = glmerControl(
                                optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.10.alt <- glmer(transfer ~ 
                                 alt_diff_belief * monitoring_ability + 
                                 role_first_half +
                                 age + 
                                 female + 
                                 christian + 
                                 muslim + 
                                 beautytask_choice + 
                                 crt_total_correct + 
                                 altruism + 
                                 positive_reciprocity + 
                                 round + 
                                 (1|subject_id),
                               family = "binomial", data = lab.data.6rounds, 
                               control = glmerControl(
                                 optimizer ='optimx', optCtrl=list(method='nlminb')))

coef.7.alt <- coeftest(reg.lab.7.alt)
coef.8.alt <- coeftest(reg.lab.8.alt)

sink("table_Appendix_E_10.tex")
stargazer(coef.7.alt, coef.8.alt, 
          reg.lab.9.alt, reg.lab.10.alt, 
          title = "Replication of Models 1 to 4 with alternative measure of a belief along a continous sacle",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Diff. in adjusted beliefs", "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism", 
                                "Positive reciprocity", "Round",
                                "Diff. in adjusted beliefs X Private transfer offered" , 
                                "Diff. in adjusted beliefs X Monitoring capacity"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(2, 2), 
          add.lines = list(c("Observations", "360", "360", "360", "360")), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 4 use random effect regressions (ologit for Models 1 and 2 and logit for Models 3 and 4); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()

response.level.labs<- c("Voting for Programmatic candidate", "Abstention", "Voting for Clientelistic candidate")
names(response.level.labs) <- c("1", "2", "3")
p1.alt<- plot_model(reg.lab.8.alt, type = "pred", terms = c("alt_diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "partner_transfer"), 
                    colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Private transfer",
                    axis.title = "Predicted probability", show.legend = TRUE) + theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1,0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") + 
  facet_grid(. ~ response.level, labeller = labeller(response.level = response.level.labs)) + theme(strip.text.x = element_text(size = 10))
p1.alt 
p1.alt[1]
ggsave("marginaleffect1alt.png", width=7, height=4, scale=1.1)

p2.alt<- plot_model(reg.lab.10.alt , type = "pred", terms = c("alt_diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "monitoring_ability"), 
                    colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Monitoring capacity",
                    axis.title = "Predicted probability", show.legend = TRUE) +  theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1, 0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") #+ geom_point() 
p2.alt
p2.alt[1]
ggsave("marginaleffect2alt.png", width=7, height=4, scale=1.1)


#####################
### Appendix E.11 ###
##################### 
reg.lab.15 <- clmm(vote_choice ~   
                            diff_belief +  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            round + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.15_2 <- clmm(vote_choice ~   
                            diff_belief *  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            round + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.16 <- glmer(transfer ~ 
                             diff_belief + monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             round + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.16_2 <- glmer(transfer ~ 
                             diff_belief * monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             round + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

coef.15 <- coeftest(reg.lab.15)
coef.15_2 <- coeftest(reg.lab.15_2)

sink("table_Appendix_E_11.tex")
stargazer(coef.15, coef.15_2, 
          reg.lab.16, reg.lab.16_2, 
          title = "Replication of Table 3 without post-manipulation variables",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Diff. in adjusted beliefs", "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Round",
                                "Diff. in adjusted beliefs X Transfer offered" , 
                                "Diff. in adjusted beliefs X Monitoring capacity"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(2, 2), 
          add.lines = list(c("Observations", "360", "360", "360", "360")), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 4 use random effect regressions (ologit for Models 1 and 2 and logit for Models 3 and 4); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()

p5<- plot_model(reg.lab.15_2, type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "partner_transfer"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Private transfer",
                axis.title = "Predicted probability", show.legend = TRUE) + theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1,0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") + 
  #geom_point() +
  facet_grid(. ~ response.level, labeller = labeller(response.level = response.level.labs)) + theme(strip.text.x = element_text(size = 10))
p5
p5[1]
ggsave("marginaleffect5.png", width=7, height=4, scale=1.1)

p6<- plot_model(reg.lab.16_2, type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "monitoring_ability"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Monitoring capacity",
                axis.title = "Predicted probability", show.legend = TRUE) +  theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1, 0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") #+ geom_point() 
p6
p6[1]
ggsave("marginaleffect6.png", width=7, height=4, scale=1.1)


####################
### Appendix E.12 ##
####################
reg.lab.17 <- clmm(vote_choice ~   
                            diff_belief +  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            crt_total_correct + 
                            altruism + 
                            positive_reciprocity + 
                            round + 
                            lagged_partner_transfer +
                            lagged_vote_choice +
                            lagged_sanctioned + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.17_2 <- clmm(vote_choice ~   
                            diff_belief *  partner_transfer + 
                            role_first_half +
                            age + 
                            female + 
                            christian + 
                            muslim + 
                            beautytask_choice + 
                            crt_total_correct + 
                            altruism + 
                            positive_reciprocity + 
                            round + 
                            lagged_partner_transfer +
                            lagged_vote_choice +
                            lagged_sanctioned + 
                            (1|subject_id),
                          data = lab.data.6rounds )

reg.lab.18 <- glmer(transfer ~ 
                             diff_belief + monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             crt_total_correct + 
                             altruism + 
                             positive_reciprocity + 
                             round + 
                             lagged_monitoring_ability +
                             lagged_transfer +
                             lagged_partner_sanctioned + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

reg.lab.18_2<- glmer(transfer ~ 
                             diff_belief * monitoring_ability + 
                             role_first_half +
                             age + 
                             female + 
                             christian + 
                             muslim + 
                             beautytask_choice + 
                             crt_total_correct + 
                             altruism + 
                             positive_reciprocity + 
                             round + 
                             lagged_monitoring_ability +
                             lagged_transfer +
                             lagged_partner_sanctioned + 
                             (1|subject_id),
                           family = "binomial", data = lab.data.6rounds, 
                           control = glmerControl(
                             optimizer ='optimx', optCtrl=list(method='nlminb')))

coef.17 <- coeftest(reg.lab.17)
coef.17_2 <- coeftest(reg.lab.17_2)

sink("table_Appendix_E_12.tex")
stargazer(coef.17, coef.17_2, 
          reg.lab.18, reg.lab.18_2, 
          title = "Replication of Table 3 with lagged variables",
          type="latex", 
          omit = c("mu_1", "sigma"), 
          covariate.labels = c( "cutoff1", "cutoff2",
                                "Diff. in adjusted beliefs", "Private transfer offered", "Monitoring capacity", 
                                "Role in the first six rounds","Age", "Female", "Christian", "Muslim", 
                                "Choice in the beauty contest", "Num. of correct answers in the CRT", "Altruism", 
                                "Positive reciprocity", "Round", 
                                "Private transfer offer in the previous round", "Vote choice in the previous round", "Punished in the previous round",
                                "Diff. in adjusted beliefs X Private transfer offered",
                                "Monitoring capacity in the previous round", "Offered a private transfer in the previous round",  "Voter was punished in the previous round", 
                                "Diff. in adjusted beliefs X Monitoring capacity"
          ),
          model.names = FALSE, 
          dep.var.labels.include = FALSE,
          column.labels   = c("Vote choice (-1/0/1)", "Transfer offer (0/1)"),
          column.separate = c(2, 2), 
          add.lines = list(c("Observations", "300", "300", "300", "300")), 
          digits = 3,
          star.cutoffs = c(0.05, 0.01), 
          omit.stat = c("rsq", "f", "aic", "ser"), 
          notes        = "Models 1 to 4 use random effect regressions (ologit for Models 1 and 2 and logit for Models 3 and 4); standard errors in parentheses; cutoffs are suppressed in the report; * p<0.05, ** p<0.01", 
          notes.append = FALSE, 
          notes.align = "l")
sink()


p7<- plot_model(reg.lab.17_2, type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "partner_transfer"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Private transfer",
                axis.title = "Predicted probability", show.legend = TRUE) + theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1,0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") + 
  facet_grid(. ~ response.level, labeller = labeller(response.level = response.level.labs)) + theme(strip.text.x = element_text(size = 10))
p7 
p7[1]
ggsave("marginaleffect7.png", width=7, height=4, scale=1.1)


p8<- plot_model(reg.lab.18_2, type = "pred", terms = c("diff_belief[-1,-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1]", "monitoring_ability"), 
                colors="bw",  ci.lvl= 0.95, title = "", legend.title = "Monitoring capacity",
                axis.title = "Predicted probability", show.legend = TRUE) +  theme_bw() + 
  coord_cartesian(expand = TRUE, ylim=c(0, 1)) + scale_x_continuous(limits = c(-1, 1), breaks = c(-1, 0, 1)) + xlab("Belief updating") + theme(legend.position="bottom") #+ geom_point() 
p8
p8[1]
ggsave("marginaleffect8.png", width=7, height=4, scale=1.1)







